clear; clc; close all;

%% 参数设置
alpha_list = linspace(2, 3, 30);
q10_list = linspace(0, 1, 30);

T = 500;
dt = 0.01;
tspan = 0:dt:T;
t_trans = T/2;

N_alpha = length(alpha_list);
N_q10 = length(q10_list);
N_total = N_alpha * N_q10;

ClassMapVec = zeros(N_total, 1);

%% 并行计算
parfor idx = 1:N_total
    ia = mod(idx-1, N_alpha) + 1;
    iq = floor((idx-1) / N_alpha) + 1;

    alpha = alpha_list(ia);
    q10 = q10_list(iq);

    % 初始条件
    x0 = [1e-6; 0; 0; q10; 0; 0];

    % 积分轨迹
    [t, X] = ode45(@(t,x) mCNN(x, alpha), tspan, x0);

    idx_steady = t > t_trans;
    x1_steady = X(idx_steady, 1);

    % 估计最大李雅普诺夫指数
    LE_max = estimate_max_LE(@(x) mCNN(x, alpha), x0, dt, T, t_trans);

    % 分类
    class_num = Classify(x1_steady, 10, LE_max);

    ClassMapVec(idx) = class_num;
end

% 重塑为二维矩阵
ClassMap = reshape(ClassMapVec, [N_q10, N_alpha]);

%% 绘图
figure;
imagesc(alpha_list, q10_list, ClassMap);
set(gca, 'YDir', 'normal');

colormap([ ...
    0.2 0.4 1.0;   % P0 稳定点 深蓝
    0.2 0.6 1.0;   % P1
    0.4 0.6 0.9;   % P2
    0.6 0.6 1.0;   % P3
    0.6 0.7 1.0;   % P4
    0.7 1.0 0.7;   % P5
    1.0 1.0 0.5;   % P6
    1.0 0.6 0.0;   % MP 复杂多周期 橙色
    0.6 0.2 0.2;   % CH 混沌 红色
    0.5 0.5 0.5    % UB 不稳定 灰色
    ]); 

caxis([0.5 10.5]);
colorbar('Ticks', 1:10, ...
    'TickLabels', {'P0','P1','P2','P3','P4','P5','P6','MP','CH','UB'});

xlabel('\alpha');
ylabel('q_{10}');
title('Fig.3 分类图');

%% 最大李雅普诺夫指数估计函数
function LE_max = estimate_max_LE(f, x0, dt, t_total, t_trans)
    delta0 = 1e-7;
    x1 = x0;
    x2 = x0 + delta0 * randn(size(x0));
    steps = round(t_total / dt);
    trans_steps = round(t_trans / dt);

    sum_ln = 0;
    count = 0;

    for i = 1:steps
        k1x1 = f(x1);
        k1x2 = f(x2);

        x1 = x1 + dt * k1x1;
        x2 = x2 + dt * k1x2;

        dist = norm(x2 - x1);
        if dist == 0
            dist = 1e-16;
        end

        if i > trans_steps
            sum_ln = sum_ln + log(dist / delta0);
            count = count + 1;
        end

        diff_vec = (x2 - x1) / dist;
        x2 = x1 + delta0 * diff_vec;
    end

    LE_max = sum_ln / (count * dt);
end

%% 系统函数
function dx = mCNN(x, alpha)
    A = [-alpha, 0, 0;
          0,     0, 0;
          0,     0, 0];
    B = [1, -4, -3.5;
         0,  1,  2;
        -1, -4, 1.5];
    X = x(1:3);
    Q = x(4:6);
    M = exp(-Q);
    Y = M .* X;
    dXdt = -X + A*Y + B*X;
    dQdt = X;
    dx = [dXdt; dQdt];
end

%% 分类函数
function class_num = Classify(traj, max_period, LE_max) 
    if max(traj) - min(traj) < 0.4
        class_num = 1; % P0 定点
        return;
    end
    [pks, ~] = findpeaks(traj);
    if length(pks) < 2
        class_num = 1; % 定点或无振荡
        return;
    end
    pks_norm = pks / max(abs(pks));
    for period = 1:max_period
        repeated = true;
        for i = 1:(length(pks_norm) - period)
            if abs(pks_norm(i) - pks_norm(i + period)) > 0.2
                repeated = false;
                break;
            end
        end
        if repeated
            class_num = period + 1; % P1为2，依次类推
            return;
        end
    end
    if LE_max > 0
        class_num = 9; % CH 混沌
    else
        class_num = 8; % MP 多周期
    end
end
